PFG 5005 - Estudo Dirigido 1: Seções de Poincaré para a Hamiltoniana de Hénon-Heiles¶
Neste notebook, implementaremos um método de integração numérica simplético para a Hamiltoniana de Hénon-Heiles e o algoritmo de Hénon para construir seções de Poincaré.
5.1. Equações do Método de Euler Simplético¶
A Hamiltoniana de Hénon-Heiles é dada por:
$$H=\frac{1}{2}(p_{1}^{2}+p_{2}^{2}+q_{1}^{2}+q_{2}^{2})+q_{1}^{2}q_{2}-\frac{1}{3}q_{2}^{3}.$$
Para aplicar o método de Euler simplético, precisamos das equações de movimento, que são as derivadas da Hamiltoniana em relação às posições e momentos.
$$\dot{q_1} = \frac{\partial H}{\partial p_1} = p_1$$ $$\dot{q_2} = \frac{\partial H}{\partial p_2} = p_2$$ $$\dot{p_1} = -\frac{\partial H}{\partial q_1} = -(q_1 + 2q_1q_2)$$ $$\dot{p_2} = -\frac{\partial H}{\partial q_2} = -(q_2 + q_1^2 - q_2^2)$$
import numpy as np
import matplotlib.pyplot as plt
import os
import warnings
# converte todos os RuntimeWarning do numpy em exceções
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(over='raise', invalid='raise', divide='raise', under='ignore')
def H(q1, q2, p1, p2):
"""
Função da Hamiltoniana de Hénon-Heiles.
"""
return 0.5 * (p1**2 + p2**2 + q1**2 + q2**2) + q1**2 * q2 - (1/3) * q2**3
5.2. Implementação do método de Euler simplético e algoritmo de Henon¶
Método de Euler¶
O método de Euler Simplético para esta Hamiltoniana separável é um mapa que atualiza as posições e momentos em duas etapas:
Atualização das Posições: $$q_1^{(n+1)}=q_1^{(n)}+\Delta t p_1^{(n)}$$ $$q_2^{(n+1)}=q_2^{(n)}+\Delta t p_2^{(n)}$$
Atualização dos Momentos (usando as novas posições): $$p_1^{(n+1)}=p_1^{(n)}-\Delta t (q_1^{(n+1)} + 2q_1^{(n+1)}q_2^{(n+1)})$$ $$p_2^{(n+1)}=p_2^{(n)}-\Delta t (q_2^{(n+1)} + (q_1^{(n+1)})^2 - (q_2^{(n+1)})^2)$$
def euler_symplectic_step(q1, q2, p1, p2, dt):
"""
Realiza um passo de integração usando o método de Euler Simplético.
"""
# Atualiza as posições
q1_new = q1 + dt * p1
q2_new = q2 + dt * p2
# Atualiza os momentos usando as novas posições
p1_new = p1 - dt * (q1_new + 2 * q1_new * q2_new)
p2_new = p2 - dt * (q2_new + q1_new**2 - q2_new**2)
return q1_new, q2_new, p1_new, p2_new
Algoritmo de Henon¶
O algoritmo de Henon encontra o ponto exato onde a trajetória cruza a seção de Poincaré $q_1 = 0$. O algoritmo de Hénon envolve um segundo passo de integração para encontrar o ponto exato de interseção com a superfície de seção. Para isso, reformulamos as equações de movimento, usando $q_1$ como a variável independente:
$$\frac{dq_2}{dq_1} = \frac{\dot{q_2}}{\dot{q_1}} = \frac{p_2}{p_1}$$ $$\frac{dp_2}{dq_1} = \frac{\dot{p_2}}{\dot{q_1}} = \frac{-(q_2 + q_1^2 - q_2^2)}{p_1}$$
Verifica o cruzamento
Checa se $q_1$ mudou de sinal (foi de negativo para positivo) e se $p_1 \ge 0$. Isso garante que a trajetória realmente passou pela seção na direção certa.Calcula o quanto falta em $q_1$
A distância até a seção é simplesmente
$$ \Delta q_1 = -q_1^{(n)} . $$Avança as outras variáveis em função de $q_1$ Em vez de usar o tempo $t$, o algoritmo usa $q_1$ como “passo”.
Retorna o ponto de interseção
O par $(q_2^{*}, p_2^{*})$ é o ponto preciso na seção $q_1 = 0$.
def henon_alg(q1_old, q2_old, p1_old, p2_old, q1, p1):
# Condição de cruzamento
if q1_old < 0 and q1 >= 0 and p1 >= 0:
# Δq1 necessário para alcançar q1=0
d_q1 = -q1_old
# Euler em função de q1
q2_intersec = q2_old + d_q1 * (p2_old / p1_old)
p2_intersec = p2_old + d_q1 * (-(q2_old + q1_old**2 - q2_old**2) / p1_old)
return q2_intersec, p2_intersec
5.3. Seção de Poincaré¶
Esta célula contém a lógica principal para a simulação. A função get_poincare_points referencia o integrador de Euler Simplético e o algoritmo de Hénon para detectar e registrar os pontos da seção de Poincaré.
A seção de Poincaré será construída no plano $(q_2, p_2)$ para a superfície de seção $q_1=0$, com a condição adicional de $p_1 \ge 0$.
def get_poincare_points(initial_state, E_target, dt, num_steps):
# (código da sua função)
q1, q2, p1, p2 = initial_state
poincare_points_q2 = []
poincare_points_p2 = []
# Loop de integração
i=0
for n in range(num_steps):
# Armazena o estado anterior para detecção de cruzeiro
q1_old, q2_old, p1_old, p2_old = q1, q2, p1, p2
# Passo de integração principal
try:
q1, q2, p1, p2 = euler_symplectic_step(q1_old, q2_old, p1_old, p2_old, dt)
except FloatingPointError:
i += 1
break # Aborta essa trajetória
# Verifica a condição de cruzamento
hit = henon_alg(q1_old, q2_old, p1_old, p2_old, q1, p1)
if hit is not None:
q2_intersec, p2_intersec = hit
poincare_points_q2.append(q2_intersec)
poincare_points_p2.append(p2_intersec)
return poincare_points_q2, poincare_points_p2, i
5.3. Seção de Poincaré para E=0.08333¶
Para construir a seção de Poincaré, precisamos de uma condição inicial que satisfaça a energia total $E=0.08333$. Vamos escolher uma condição inicial com $q_1=0$ para começar, e então usar a equação da Hamiltoniana para encontrar o valor de $p_2$.
$H = E \implies \frac{1}{2}(p_{1}^{2}+p_{2}^{2}+q_{1}^{2}+q_{2}^{2})+q_{1}^{2}q_{2}-\frac{1}{3}q_{2}^{3} = E$
Fixando $q_1=0$, $p_1=0.2$ e $q_2=0$, podemos resolver para $p_2$.
$E = 0.08333$
$0.08333 = \frac{1}{2}(0.2^2 + p_2^2 + 0^2 + 0^2) + 0 - 0$ $0.08333 = \frac{1}{2}(0.04 + p_2^2)$ $0.16666 = 0.04 + p_2^2$ $p_2^2 = 0.12666 \implies p_2 \approx \pm 0.35589$
Vamos usar $p_2 = 0.35589$.
Figura 4 da Ref. [4]
# Parâmetros da simulação
E_target = 0.08333
dt = 0.01
num_steps = 10**6
# Geração de múltiplas condições iniciais em uma grade
# Aumente o número de pontos para q2_vals e p2_vals para ter mais trajetórias
# Por exemplo, uma grade 30x30 resulta em 900 combinações
q2_vals = np.linspace(-0.5, 0.5, 30)
p2_vals = np.linspace(-0.5, 0.5, 30)
total_trajectories = len(q2_vals) * len(p2_vals)
all_q2_points = []
all_p2_points = []
num_valid_initials = 0
trajectory_counter = 0
print(f"Iniciando simulações para E = {E_target}...")
# O loop principal para gerar as múltiplas trajetórias
for q2_0 in q2_vals:
for p2_0 in p2_vals:
trajectory_counter += 1
# A condição inicial de p1 é calculada a partir da energia total
# H = 1/2(p1^2 + p2^2 + q1^2 + q2^2) + q1^2*q2 - 1/3*q2^3
# Para q1=0, a equação da energia se torna E = 1/2(p1^2 + p2^2 + q2^2) - 1/3*q2^3
# E pode ser reescrita para isolar p1: p1^2 = 2E - (p2^2 + q2^2) + (2/3)q2^3
p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3
if p1_sq >= 0:
p1_calc = np.sqrt(p1_sq)
initial_state = (0.0, q2_0, p1_calc, p2_0)
# print(f" > E = {E_target}: Trajetória {trajectory_counter}/{total_trajectories}...")
q2_points, p2_points = get_poincare_points(initial_state, E_target, dt, num_steps)
all_q2_points.extend(q2_points)
all_p2_points.extend(p2_points)
num_valid_initials += 1
# Bloco para plotar o gráfico
if num_valid_initials > 0:
plt.figure(figsize=(10, 10))
plt.scatter(all_q2_points, all_p2_points, s=0.01, alpha=0.5)
# Define os limites de visualização fixos
q2_min, q2_max = -0.6, 0.6
p2_min, p2_max = -0.6, 0.6
plt.xlim(q2_min, q2_max)
plt.ylim(p2_min, p2_max)
plt.xlabel('$q_2$', fontsize=12)
plt.ylabel('$p_2$', fontsize=12)
plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
plt.grid(True)
plt.show()
print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n")
Iniciando simulações para E = 0.08333...
Concluído para E = 0.08333! Total de pontos: 718199
# Parâmetros da simulação
E_target = 0.125
dt = 0.01
num_steps = 10**6
# Geração de múltiplas condições iniciais em uma grade
# Aumente o número de pontos para q2_vals e p2_vals para ter mais trajetórias
# Por exemplo, uma grade 30x30 resulta em 900 combinações
q2_vals = np.linspace(-0.5, 0.5, 30)
p2_vals = np.linspace(-0.5, 0.5, 30)
total_trajectories = len(q2_vals) * len(p2_vals)
all_q2_points = []
all_p2_points = []
num_valid_initials = 0
trajectory_counter = 0
print(f"Iniciando simulações para E = {E_target}...")
# O loop principal para gerar as múltiplas trajetórias
for q2_0 in q2_vals:
for p2_0 in p2_vals:
trajectory_counter += 1
# A condição inicial de p1 é calculada a partir da energia total
# H = 1/2(p1^2 + p2^2 + q1^2 + q2^2) + q1^2*q2 - 1/3*q2^3
# Para q1=0, a equação da energia se torna E = 1/2(p1^2 + p2^2 + q2^2) - 1/3*q2^3
# E pode ser reescrita para isolar p1: p1^2 = 2E - (p2^2 + q2^2) + (2/3)q2^3
p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3
if p1_sq >= 0:
p1_calc = np.sqrt(p1_sq)
initial_state = (0.0, q2_0, p1_calc, p2_0)
# print(f" > E = {E_target}: Trajetória {trajectory_counter}/{total_trajectories}...")
q2_points, p2_points = get_poincare_points(initial_state, E_target, dt, num_steps)
all_q2_points.extend(q2_points)
all_p2_points.extend(p2_points)
num_valid_initials += 1
# Bloco para plotar o gráfico
if num_valid_initials > 0:
plt.figure(figsize=(10, 10))
plt.scatter(all_q2_points, all_p2_points, s=0.01, alpha=0.5)
# Define os limites de visualização fixos
q2_min, q2_max = -0.6, 0.7
p2_min, p2_max = -0.6, 0.7
plt.xlim(q2_min, q2_max)
plt.ylim(p2_min, p2_max)
plt.xlabel('$q_2$', fontsize=12)
plt.ylabel('$p_2$', fontsize=12)
plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
plt.grid(True)
plt.show()
print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n")
Iniciando simulações para E = 0.125...
Concluído para E = 0.125! Total de pontos: 1009720
5.5. Séries temporais¶
Séries temporais de $q_1$, $p_1$, $q_2$, $q_2$, para E=0.125 de trajetórias periódicas.
import numpy as np
import matplotlib.pyplot as plt
# Hamiltoniana Hénon–Heiles (para monitorar energia, opcional)
def H(q1, q2, p1, p2):
return 0.5*(p1*p1 + p2*p2 + q1*q1 + q2*q2) + q1*q1*q2 - (1.0/3.0)*q2**3
# Euler simplético (Euler–Cromer): atualiza q com p^n e p com grad V em q^{n+1}
def euler_symplectic_step(q1, q2, p1, p2, dt):
# q^{n+1}
q1_new = q1 + dt*p1
q2_new = q2 + dt*p2
# grad V(q1,q2) = ( q1 + 2*q1*q2 , q2 + q1^2 - q2^2 )
p1_new = p1 - dt*(q1_new + 2.0*q1_new*q2_new)
p2_new = p2 - dt*(q2_new + q1_new**2 - q2_new**2)
return q1_new, q2_new, p1_new, p2_new
# gera séries temporais
def simulate_time_series(initial_state, dt, num_steps):
q1, q2, p1, p2 = initial_state
t = np.arange(num_steps+1)*dt
q1_series = np.empty(num_steps+1)
q2_series = np.empty(num_steps+1)
p1_series = np.empty(num_steps+1)
p2_series = np.empty(num_steps+1)
q1_series[0], q2_series[0], p1_series[0], p2_series[0] = q1, q2, p1, p2
for n in range(num_steps):
q1, q2, p1, p2 = euler_symplectic_step(q1, q2, p1, p2, dt)
q1_series[n+1], q2_series[n+1], p1_series[n+1], p2_series[n+1] = q1, q2, p1, p2
return t, q1_series, p1_series, q2_series, p2_series
def plot_time_series(t, q1, p1, q2, p2, E_label=None):
fig, axs = plt.subplots(2, 2, figsize=(11, 6), sharex=True)
axs = axs.ravel()
axs[0].plot(t, q1, lw=1); axs[0].set_ylabel(r"$q_1$")
axs[1].plot(t, p1, lw=1); axs[1].set_ylabel(r"$p_1$")
axs[2].plot(t, q2, lw=1); axs[2].set_ylabel(r"$q_2$")
axs[3].plot(t, p2, lw=1); axs[3].set_ylabel(r"$p_2$")
axs[3].set_xlabel("t")
axs[2].set_xlabel("t")
title = "Séries temporais (Euler simplético)"
if E_label is not None: title += f" — E = {E_label}"
fig.suptitle(title, y=0.98, fontsize=14)
for ax in axs: ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()
E_1 = 0.125
dt = 1e-3 # menor para séries suaves
num_steps = 200000 # ajuste para cobrir tempo suficiente
q1_0, q2_0 = 0.0, 0.0
p2_0 = np.sqrt(0.21) # ~0.45826
# p1 a partir da energia (para q1=0)
p1_0 = np.sqrt(2*E_1 - (p2_0**2 + q2_0**2) + (2/3.0)*q2_0**3)
initial_state = (q1_0, q2_0, p1_0, p2_0)
t, q1_t, p1_t, q2_t, p2_t = simulate_time_series(initial_state, dt, num_steps)
plot_time_series(t, q1_t, p1_t, q2_t, p2_t, E_label=E_1)
print("Energia inicial:", H(*initial_state))
print("Energia final :", H(q1_t[-1], q2_t[-1], p1_t[-1], p2_t[-1]))
Energia inicial: 0.125 Energia final : 0.12504638709312946
5.6. Séries temporais¶
Séries temporais para E=0.125, trajetória caótica.
import numpy as np
import matplotlib.pyplot as plt
def H(q1, q2, p1, p2):
return 0.5*(p1*p1 + p2*p2 + q1*q1 + q2*q2) + q1*q1*q2 - (1.0/3.0)*q2**3
def euler_symplectic_step(q1, q2, p1, p2, dt):
q1_new = q1 + dt*p1
q2_new = q2 + dt*p2
# grad V(q) = ( q1 + 2*q1*q2 , q2 + q1**2 - q2**2 )
p1_new = p1 - dt*(q1_new + 2.0*q1_new*q2_new)
p2_new = p2 - dt*(q2_new + q1_new**2 - q2_new**2)
return q1_new, q2_new, p1_new, p2_new
def simulate_time_series(initial_state, dt, num_steps):
q1, q2, p1, p2 = initial_state
t = np.arange(num_steps+1)*dt
q1s = np.empty(num_steps+1); q2s = np.empty(num_steps+1)
p1s = np.empty(num_steps+1); p2s = np.empty(num_steps+1)
q1s[0], q2s[0], p1s[0], p2s[0] = q1, q2, p1, p2
for n in range(num_steps):
q1, q2, p1, p2 = euler_symplectic_step(q1, q2, p1, p2, dt)
q1s[n+1], q2s[n+1], p1s[n+1], p2s[n+1] = q1, q2, p1, p2
return t, q1s, p1s, q2s, p2s
E = 0.16667
q1_0, q2_0 = 0.0, 0.0
p1_0 = 0.2
p2_0 = np.sqrt(0.29334) # ≈ 0.5416 (compatível com E=0.16667)
dt = 1e-3
num_steps = 300_000
t, q1_t, p1_t, q2_t, p2_t = simulate_time_series((q1_0, q2_0, p1_0, p2_0), dt, num_steps)
fig, axs = plt.subplots(2,2, figsize=(11,6), sharex=True)
axs = axs.ravel()
axs[0].plot(t, q1_t, lw=0.8); axs[0].set_ylabel(r"$q_1$")
axs[1].plot(t, p1_t, lw=0.8); axs[1].set_ylabel(r"$p_1$")
axs[2].plot(t, q2_t, lw=0.8); axs[2].set_ylabel(r"$q_2$"); axs[2].set_xlabel("t")
axs[3].plot(t, p2_t, lw=0.8); axs[3].set_ylabel(r"$p_2$"); axs[3].set_xlabel("t")
fig.suptitle(f"Séries temporais — trajetória caótica (E={E})", y=0.98)
for ax in axs: ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()
E_2 = 0.125
dt = 0.01
num_steps = 10**6
# Condição inicial para E = 0.16667
# 0.16667 = 0.5 * (0.2**2 + p_2**2)
# 0.33334 = 0.04 + p_2**2
# p_2**2 = 0.29334 => p_2 = sqrt(0.29334) ~ 0.5416
q1_0, q2_0, p1_0, p2_0 = 0.0, 0.0, 0.2, np.sqrt(0.29334)
initial_state = (q1_0, q2_0, p1_0, p2_0)
q2_points, p2_points = get_poincare_points(initial_state, E_2, dt, num_steps)
# Plotar os resultados
plt.figure(figsize=(8, 8))
plt.scatter(q2_points, p2_points, s=1)
plt.xlabel('$q_2$', fontsize=12)
plt.ylabel('$p_2$', fontsize=12)
plt.title(f'Seção de Poincaré (E = {E_2})', fontsize=14)
plt.grid(True)
plt.show()
print(f"Número de pontos na seção: {len(q2_points)}")
Número de pontos na seção: 1414
5.6. Seção de Poincaré para E=0.16667¶
# Parâmetros da simulação
E_target = 0.16667
dt = 0.01
num_steps = 10**6
# Geração de múltiplas condições iniciais em uma grade
# Aumente o número de pontos para q2_vals e p2_vals para ter mais trajetórias
# Por exemplo, uma grade 30x30 resulta em 900 combinações
q2_vals = np.linspace(-0.5, 0.5, 30)
p2_vals = np.linspace(-0.5, 0.5, 30)
total_trajectories = len(q2_vals) * len(p2_vals)
all_q2_points = []
all_p2_points = []
num_valid_initials = 0
trajectory_counter = 0
print(f"Iniciando simulações para E = {E_target}...")
# O loop principal para gerar as múltiplas trajetórias
i=0
for q2_0 in q2_vals:
for p2_0 in p2_vals:
trajectory_counter += 1
# A condição inicial de p1 é calculada a partir da energia total
# H = 1/2(p1^2 + p2^2 + q1^2 + q2^2) + q1^2*q2 - 1/3*q2^3
# Para q1=0, a equação da energia se torna E = 1/2(p1^2 + p2^2 + q2^2) - 1/3*q2^3
# E pode ser reescrita para isolar p1: p1^2 = 2E - (p2^2 + q2^2) + (2/3)q2^3
p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3
if p1_sq >= 0:
p1_calc = np.sqrt(p1_sq)
initial_state = (0.0, q2_0, p1_calc, p2_0)
# print(f" > E = {E_target}: Trajetória {trajectory_counter}/{total_trajectories}...")
q2_points, p2_points, j = get_poincare_points(initial_state, E_target, dt, num_steps)
i += j
all_q2_points.extend(q2_points)
all_p2_points.extend(p2_points)
num_valid_initials += 1
# Bloco para plotar o gráfico
if num_valid_initials > 0:
plt.figure(figsize=(10, 10))
plt.scatter(all_q2_points, all_p2_points, s=0.01, alpha=0.5)
# Define os limites de visualização fixos
q2_min, q2_max = -0.6, 1.0
p2_min, p2_max = -0.6, 1.0
plt.xlim(q2_min, q2_max)
plt.ylim(p2_min, p2_max)
plt.xlabel('$q_2$', fontsize=12)
plt.ylabel('$p_2$', fontsize=12)
plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
plt.grid(True)
plt.show()
print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n Total de Trajetórias Abortadas: {i}\n")
Iniciando simulações para E = 0.16667...
/var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:10: RuntimeWarning: overflow encountered in scalar multiply p1_new = p1 - dt * (q1_new + 2 * q1_new * q2_new) /var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:11: RuntimeWarning: overflow encountered in scalar power p2_new = p2 - dt * (q2_new + q1_new**2 - q2_new**2) /var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:11: RuntimeWarning: invalid value encountered in scalar subtract p2_new = p2 - dt * (q2_new + q1_new**2 - q2_new**2) /var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:10: RuntimeWarning: invalid value encountered in scalar subtract p1_new = p1 - dt * (q1_new + 2 * q1_new * q2_new)
Concluído para E = 0.16667! Total de pontos: 979611 Total de Trajetórias Abortadas: 0
Seção de Poincaré para E variando de 0.0050 a 0.1667¶
Vamos fazer a mesma coisa para cada energia na lista: 0.0050, 0.0070, 0.0089, 0.0100, 0.0117, 0.0133, 0.0148, 0.0167, 0.0183, 0.0188, 0.0217, 0.0250, 0.0267, 0.0300, 0.0333, 0.0350, 0.0383, 0.0417, 0.0433, 0.0467, 0.0500, 0.0517, 0.0550, 0.0583, 0.0600, 0.0633, 0.0667, 0.0683, 0.0717, 0.0750, 0.0767, 0.0800, 0.0833, 0.0850, 0.0883, 0.0917, 0.0933, 0.0967, 0.1000, 0.1017, 0.1050, 0.1083, 0.1100, 0.1133, 0.1167, 0.1183, 0.1217, 0.1250, 0.1267, 0.1300, 0.1333, 0.1350, 0.1383, 0.1417, 0.1433, 0.1467, 0.1500, 0.1517, 0.1550, 0.1583, 0.1600, 0.1633, 0.1667 , variando também as trajetórias em cada energia para 50 valores de q e 50 valores de p gerando 2500 condições iniciais distintas.
E_list = [
0.0050, 0.0070, 0.0089, 0.0100, 0.0117, 0.0133, 0.0148, 0.0167, 0.0183,
0.0188, 0.0217, 0.0250, 0.0267, 0.0300, 0.0333, 0.0350, 0.0383, 0.0417,
0.0433, 0.0467, 0.0500, 0.0517, 0.0550, 0.0583, 0.0600, 0.0633, 0.0667,
0.0683, 0.0717, 0.0750, 0.0767, 0.0800, 0.0833, 0.0850, 0.0883, 0.0917,
0.0933, 0.0967, 0.1000, 0.1017, 0.1050, 0.1083, 0.1100, 0.1133, 0.1167,
0.1183, 0.1217, 0.1250, 0.1267, 0.1300, 0.1333, 0.1350, 0.1383, 0.1417,
0.1433, 0.1467, 0.1500, 0.1517, 0.1550, 0.1583, 0.1600, 0.1633, 0.1667
]
# Parâmetros da simulação
dt = 0.01
num_steps = 10**6
# Geração de múltiplas condições iniciais em uma grade
# q2_vals = np.linspace(-0.5, 0.5, 50)
# p2_vals = np.linspace(-0.5, 0.5, 50)
# total_trajectories = len(q2_vals) * len(p2_vals)
# Geração de múltiplas condições iniciais em uma grade
q2_vals_grid = np.linspace(-0.5, 0.5, 50)
p2_vals_grid = np.linspace(-0.5, 0.5, 50)
q2_p2_combinations = np.array(np.meshgrid(q2_vals_grid, p2_vals_grid)).T.reshape(-1, 2)
total_combinations = len(q2_p2_combinations)
# Define o número inicial e final de trajetórias
min_trajectories = 350
max_trajectories = 1000
# Cria a pasta 'imagens' para guardar os gráficos de cada seção de Poincaré
output_dir = "imagens"
if not os.path.exists(output_dir):
os.makedirs(output_dir)
print(f"Pasta '{output_dir}' criada.")
# Define limites fixos para os eixos X e Y para todos os gráficos
# Isso garante que a escala não mude, o que é ideal para um GIF.
q2_min, q2_max = -0.6, 0.6
p2_min, p2_max = -0.6, 0.6
Pasta 'imagens' criada.
# A lista de combinações selecionadas, que irá aumentar conforme aumenta a energia
selected_combinations = np.empty((0, 2))
# O conjunto de índices já presentes na simulação
simulated_indices = set()
# Itera sobre a lista de energias com um índice
for i, E_target in enumerate(E_list):
all_q2_points = []
all_p2_points = []
num_valid_initials = 0
# Calcula o número total de trajetórias que deveriam ser testadas para essa energia
num_trajectories_target = int(min_trajectories + (max_trajectories - min_trajectories) * (i / (len(E_list) - 1)))
# Adiciona novas trajetórias se a lista atual for menor que o alvo
trajectories_to_add = num_trajectories_target - len(selected_combinations)
if trajectories_to_add > 0:
# Pega os índices que ainda não foram usados
available_indices = np.array(list(set(range(total_combinations)) - simulated_indices))
if len(available_indices) > 0:
# Seleciona aleatoriamente novos índices
new_indices = np.random.choice(available_indices, min(trajectories_to_add, len(available_indices)), replace=False)
# Adiciona as novas combinações à lista
selected_combinations = np.vstack([selected_combinations, q2_p2_combinations[new_indices]])
# Adiciona os novos índices ao conjunto de usados
simulated_indices.update(new_indices)
print(f"Iniciando simulações para E = {E_target} com {len(selected_combinations)} trajetórias...")
trajectory_counter = 0
for q2_0, p2_0 in selected_combinations:
trajectory_counter += 1
p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3
if p1_sq >= 0:
p1_calc = np.sqrt(p1_sq)
initial_state = (0.0, q2_0, p1_calc, p2_0)
# print(f" > E = {E_target}: Trajetória {trajectory_counter}/{len(selected_combinations)}...")
q2_points, p2_points = get_poincare_points(initial_state, E_target, dt, num_steps)
all_q2_points.extend(q2_points)
all_p2_points.extend(p2_points)
num_valid_initials += 1
if num_valid_initials > 0:
plt.figure(figsize=(10, 10))
plt.scatter(all_q2_points, all_p2_points, s=0.02, alpha=0.5)
# Define os limites de visualização fixos
plt.xlim(q2_min, q2_max)
plt.ylim(p2_min, p2_max)
plt.xlabel('$q_2$', fontsize=12)
plt.ylabel('$p_2$', fontsize=12)
plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
plt.grid(True)
filename = f'simulacao_{E_target:.4f}'.replace('.', '_')
filepath = os.path.join(output_dir, filename)
plt.savefig(filepath, bbox_inches='tight')
print(f"Gráfico salvo como '{filepath}'")
plt.close() # Fecha a figura para não sobrecarregar a memória
# print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n")
else:
print(f"Nenhuma trajetória válida encontrada para E = {E_target}. Aumente o intervalo da grade ou a energia.\n")
Gif para as simulações¶
import glob
from PIL import Image as PILImage
from IPython.display import Image as DispImage, display
# 1) Gerar o GIF (se já tiver o .gif, pule esta parte)
frames = [Image.open(p) for p in sorted(glob.glob("imagens/*.png"))]
frames[0].save(
"animacao.gif",
save_all=True,
append_images=frames[1:],
duration=250, # ms por fram
loop=0, # 0 = loop infinito
disposal=2 # ajuda a evitar “fantasmas” de frames
)
# 2) Exibir EMBUTINDO no output (vai em base64 no HTML exportado)
# display(IPImage(filename="animacao.gif", format="gif", embed=True))
display(DispImage(filename="animacao.gif", format="gif", embed=True))
<IPython.core.display.Image object>
import base64
from IPython.display import HTML
from pathlib import Path
# Lê e converte para base64
data = Path("animacao.gif").read_bytes()
b64 = base64.b64encode(data).decode("ascii")
# Embute direto no output (vai sair no HTML final)
HTML(f'<img src="data:image/gif;base64,{b64}">')